Spontaneous thermal runaway as an ultimate failure mechanism of materials 
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The first theoretical estimate of the shear strength of a perfect crystal was given by Frenkel [Z. 
Phys. 37, 572 (1926)]. He assumed that as slip occurred, two rigid atomic rows in the crystal would 
move over each other along a slip plane. Based on this simple model, Frenkel derived the ultimate 
shear strength to be about one tenth of the shear modulus. Here we present a theoretical study 
showing that catastrophic material failure may occur below Frenkel's ultimate limit as a result of 
thermal runaway. We demonstrate that the condition for thermal runaway to occur is controlled by 
only two dimensionless variables and, based on the thermal runaway failure mechanism, we calculate 
the maximum shear strength CTc of viscoelastic materials. Moreover, during the thermal runaway 
process, the magnitude of strain and temperature progressively localize in space producing a narrow 
region of highly deformed material, i.e. a shear band. We then demonstrate the relevance of this 
new concept for material failure known to occur at scales ranging from nanometers to kilometers. 
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It is well known that the shear strength of real crys- 
tals is typically several orders of magnitude smaller than 
Frenkel's ultimate shear strength limit. This discrepancy 
is explained by the fact that real crystals contain de- 
fects such as dislocations which lower the shear strength 
dramatically. Nevertheless, some materials, like rocks 
in the Earth's interior and metallic glasses, apparently 
have strengths approaching Frenkel's theoretical limit. 
For these materials it is reasoned that the mobility of 
the defects is in one way or another reduced. For in- 
stance, the closure of cracks at high confining pressures, 
non-planar crystal structure of minerals and disorder of 
mineral grain orientations are all factors contributing to 
the high strength of rocks in the Earth's interior. The 
high strength of metallic glasses is attributed to the high 
degree of structural disorder causing dislocations to ex- 
perience a large number of obstacles, reducing their mo- 
bility and inhibiting plastic flow. 

However, even if dislocation mobility is greatly re- 
duced, materials subjected to increasing loads do not 
necessarily fail according to Frenkel's model. When sub- 
jected to a high shear stress which approaches, but is 
reproducibly lower than Frenkel's limit, these materials 
may fail by deformation localized on a single or a few 
regions (shear bands) having thicknesses that are orders 
of magnitude larger than interatomic spacing, but which 
are still very narrow compared to the deforming sam- 
ple size. Moreover, this extreme localization of shear is 
often accompanied by extensive melting and resolidifica- 
tion of the material. This mode of failure is manifested 
by the so-called pseudotachylytes in geological outcrops 
(cm thick deformation bands believed to be the traces 
of large paleoearthquakes) [l|, 0] and catastrophic failure 
along ca. 10 nm thick shear bands reported in metallic 
glasses 

The occurrence of localized shear zones in the deeper 
parts of the earth's lithosphere poses a problem due to its 
paradoxical nature. Earthquakes are usually attributed 
to the intrinsic stick-slip character of brittle failure or 
frictional sliding on a pre-existing fault. Brittle failure is 



expected as the primary mode of rock failure down to 15- 
20 kms, the so-called seismogenic zone. Due to the high 
confining pressure at depth, however, these mechanisms 
alone are not plausible explanations for earthquakes oc- 
curring significantly deeper than the seismogenic zone. 
Several potential rock-weakening mechanisms have in- 
stead been proposed to facilitate deep earthquakes, in- 
cluding structural changes such as dehydration embrit- 
tlement and olivine spinel phase transition m| or ther- 
mal softening leading to shear instability 0, HI . Recent 
assessment of these models using indirect seismological 
evidence argues in favor of a temperature-activated phe- 
nomenon, such as thermal shear instabilities, due to ap- 
parent temperature dependence of deep earthquakes [9|]. 
Similarly, structural strain softening and thermal soften- 
ing are possible explanations for shear band formation in 
metallic glasses Ifl, HI- A recent experiment how- 
ever, revealed that the heated zone associated with an 
individual band was much wider than the shear band 
itself. It was concluded [3| that shear-band operation 
could not be fully adiabatic and temperature rise could 
not therefore be the factor controlling shear-band thick- 
ness, thus rejecting the thermal softening mechanism. In 
contrast, our work presented below shows that the non- 
adiabaticity of the thermal softening process may in fact 
be the cause of strain-localization inside the hot zone, 
and we thus provide an alternative interpretation to the 
experimental results in ref. 

It is well accepted that even below the conventional 
elastic limit, most real materials show non-elastic rheo- 
logical responses such as creep under constant load and 
stress relaxation under constant extension induced by 
thermally excited defects and imperfections. In accor- 
dance with these properties, most materials may be char- 
acterized as viscoelastic, i.e. the rheology contains both 
viscous and elastic components. Since the phenomena of 
creep and relaxation are thermally activated processes, 
the viscosity is strongly temperature dependent (e.g. Ar- 
rhenius) and it is, in general, a non-linear function of the 
shear stress [12]. The strong temperature dependence of 



2 




'bg 



To h\ 



'bg 



FIG. 1: Initial setup of the viscoelastic slab model discussed in 
the text (cross-section in the a;?/- plane). The slab is in a state 
of stress of simple shear with zero velocity {v — 0) bound- 
ary conditions. The shear stress a, constant throughout the 
slab, initially (t = 0) equals the maximum value ao and sub- 
sequently decreases with time due to relaxation and viscous 
deformation in the interior. The blue and red lines show the 
initial velocity and temperature profiles, respectively. The 
shaded region illustrates a small perturbation in temperature 
To of width h at the slab center. Elsewhere, the background 
temperature is Tjg. The geometry of the strain rate profile 
concurs with that of the temperature profile. 



the viscosity has important imphcations as it leads to 
thermal softening of the material. Indeed, as first noted 
by Griggs and Baker Q, such a physical system is in- 
herently unstable: an increase in strain rate in a weaker 
zone causes a local temperature rise due to viscous dis- 
sipation and weakens the zone even further. At high 
stresses, viscous dissipation becomes substantial, and if 
heat is generated faster than it is conducted away, the 
local increase in temperature and strain rate is strongly 
amplified. Under those conditions a positive feedback 
between temperature rise and viscous dissipation is es- 
tablished and a thermal runaway develops. To determine 
whether the thermal runaway mechanism can explain the 
aforementioned material failure, we approach the prob- 
lem by considering a simple viscoelastic model which ac- 
counts for the non-elastic behavior below the ultimate 
yield point (Frenkel's limit). The temperature T in the 
model is determined by the equation for energy conser- 
vation which is coupled, through temperature dependent 
viscosity, to the rhcology equation. 

Our one-dimensional model (see fig. [1]) consists of a 
viscoelastic slab of width L at initial temperature Tbg 
except in the small central region having width h and 
slightly elevated temperature Tq. The boundaries are 
maintained at the temperature T^g. Our objective is to 
search for spontaneous modes of internal failure not aided 
or triggered by the effect of additional far-field deforma- 
tion. Hence, for time t > we impose zero velocity v 
at the boundaries and assume that, without addressing 



the loading history {t < 0), the slab initially {t = 0) is 
subjected to a shear stress ctq- The shear stress a in the 
slab satisfies the equation for conservation of momentum 



dx 



0, 



(1) 



which shows that a is independent of x and hence only 
a function of the time t. The viscoelastic rheology is 
represented by the Maxwell model [l^l, and is given by 
the equation 



dv 1 1 da 



(2) 



where v{x,t) is the velocity, G is the constant shear mod- 
ulus and /^(T, a) is the viscosity. The dependence of fi 
on T and a may be written as 



(3) 



where A and n are constants, E is the activation energy 
and R — 8.3 JK^^mole""'^ is the universal gas constant. 
Since a{t) is independent of x, it follows from eq. ([2]) 
that the geometry of the strain rate [dv/dx) profile at 
any instant concurs with that of the temperature profile 
T{x,t). Utilizing the zero velocity boundary condition, 
equation ^ may be integrated to obtain the equation 
which governs the time-dependence of a: 
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The temperature is determined by the equation for 
ergy conservation 
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where the last term accounts for viscous dissipation in the 
system. Equations (jH) and ([5]) constitute a closed system 
of coupled ordinary and partial differential equations for 
two unknown functions a(t) and T{x,t). 

First, to determine the conditions necessary for ther- 
mal runaway to occur, a linear stability analysis was car- 
ried out. Equation Q was approximated for the initial 
stages by substituting the initial conditions for T and 
a. This yields a characteristic time = /io/(2GAp) 
for stress relaxation. Here /io = /i(To,(To) and Ap — 
h/L + e^/^^o~^/^^''9 is a factor which characterizes the 
initial perturbation. Linearization of the temperature 
equation with account of stress relaxation yields that 
the growth of the perturbation in the initial stages is 
controlled by the two dimensionless variables T^/Td and 
c!o/ac, where the thermal diffusion time = j k {k is 
the thermal diffusivity) and the newly introduced stress 



aa^ \ l 2Ap ^ Tq . 



(6) 



3 




FIG. 2: Dependence of the maximum temperature rise ISTmax on the dimensionless variables oa/ac and Tr/r^. a, Contour 
plot of ATmax scaled by the adiabatic temperature rise AT^ax as a function of o"o/o"c and Tr/rd. The dark lines are contour 
lines. Due to the computational effort of fully resolving the entire self-localizing thermal runaway processes, the temperature 
at the very late stages of these processes are not presented in this plot, b, Profiles of the temperature rise AT = T — To 
inside the initially perturbed zone |2;| < h/2 at the time when the maximum temperature is reached. The solid line illustrates 
the self-localizing thermal runaway process at the location of the cross in a. The dashed line illustrates the adiabatic thermal 
runaway at the location of the dot in a. 



Here C denotes the heat capacity per volume. The so- 
lution to the linearized temperature equation is found to 
be unstable if o'o/o'c > f {^r I Td) ■ In the limit r^/rd << 1, 
which corresponds to near adiabatic conditions, the func- 
tion / quickly approaches the lower bound / ~ 1 and the 
solution therefore becomes unstable if ctq > <Jc- Thus 
CTc is the critical stress above which a thermal runaway 
may occur and therefore provides an estimate of the max- 
imum shear strength of viscoelastic materials. Initial 
stages of thermal runaway instability in a more general 
setup was recently investigated in ref . IJ] , including two- 
dimensional verification of the one dimensional predic- 
tions. The results of our linear stability analysis are in 
agreement with these numerical estimates in the limit of 
vanishing boundary velocity. 



Non-linear evolutions of the unstable runaway modes 
rapidly deviate from the exponential growth in time pre- 
dicted by linear analysis. Since important information 
about the deformation process can be inferred from the 
increase in temperature, we choose the maximum tem- 
perature rise ISTmax = Tmax — Tq during the deformation 
process as our main physical quantity to study. This en- 
ables us to quantify even the later stages of the thermal 
runaway process not considered in the linear analysis. A 
simple estimate of AT^ax during thermal runaway may 
be obtained assuming adiabatic conditions. In this case 
all the elastic energy in the system is uniformly dissipated 
as heat in the perturbed zone and overall energy balance 



yields the adiabatic temperature rise 

AT°- = 



2hGC 



(7) 



Guided by these analytical estimates, the complete time 
evolution of T and a was subsequently investigated by 
numerical methods. We simplified the problem by di- 
mensional analysis reducing it from one containing thir- 
teen dimensional parameters to one containing six di- 
mensionless parameters. The dimensionless form of the 
coupled set of equations ([4]) and ^ were solved numer- 
ically using a finite-difference method with non-uniform 
mesh and a tailored variable time-step in order to resolve 
the highly non-linear effect of localization. We have sys- 
tematically varied all six dimensionless parameters and 
computed AT„iax for each temperature evolution. Re- 
markably, it is possible to present AT^ax normalized by 
the adiabatic temperature rise (eq. ([71)) as a function 
of only two combinations of parameters, namely (Tq/cTc 
and Tr/Td, as previously suggested by the linear stabil- 
ity analysis. A representative set of runs is shown in 
fig. [Da,. This "phase diagram" was computed by vary- 
ing two of the dimensionless parameters and fixing the 
remaining four. The plot exhibits a low-temperature re- 
gion corresponding to stable deformation processes, and 
a high-temperature region corresponding to thermal run- 
away processes. These regions are sharply distinguished 
by a critical boundary having a location that correlates 
well with stability-predictions of the linear analysis. The 
phase diagram is "representative" in that it is insensitive 
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to which two out of the six dimensionless parameters are 
varied, keeping the remaining four fixed. 

In the neighborhood of the critical boundary in the 
high-temperature region, however, the temperatures are 
found to be much larger than the adiabatic temperature 
rise AT^^^ . In this region we observe a continuous local- 
ization of the temperature and strain profiles during the 
deformation process, i.e. the runaway is spatially self- 
localizing. This localization effect is illustrated for the 
temperature profile in fig. 2b. The elastic energy is thus 
dissipated in a zone much narrower than the width of 
the initial perturbation resulting in much larger tempera- 
tures. The self-localization of the runaway process arises 
from the effects of thermal diffusion: by diffusion the 
temperature profile acquires a peak in the center where 
the effect of the positive feedback mechanism accordingly 
is maximized. The runaway therefore accelerates faster 
in the center than in the regions outside and the deforma- 
tion process finally terminates in a highly localized shear 
band with a characteristic width much smaller than the 
characteristic width h of the initial perturbation. 

To evaluate the relevance of thermal runaway as a po- 
tential failure mechanism in nature, we now consider 
two case examples comparing critical stress for ther- 
mal runaway with Frenkel's ultimate shear strength limit 
(cr ~ G/IO). First, we estimate the condition necessary 
for thermal runaway to occur in olivine-dominated man- 
tle rocks. For olivine (see ref. [H) G = 7 x 10^° Pa, 
C = 3 X 10^ Jm-^K-i and E = 5,23 x 10^ Jmole-\ 
while we assume Tbg = 700 K (corresponding to a depth 
of about 40 km). To = 701-720 K and /i = 10 m, i = 10 
km. Substituting these values in equation ^ we pre- 
dict that thermal runaway possibly occurs if the stress 
exceeds a critical value in the range CTc = 0.5-1.7 GPa, 
i.e. (Tc = G/140-G/40. This estimate agrees rather well 
with the values observed in experimental studies of rock 
deformation under high confining pressure, where typical 
failure stresses are in the range 0.5-2 GPa 0,[i3. Sec- 
ond, for a metallic glass, typical values are G = 2, 6 x 10^ 
Jm-3K-i (ref. 0]), G = 34 GPa(ref. [H) and E = 
100-400 kJmole-i (ref. and Assuming h ^ 1 

/im, L — I cm, Thg — 300 K (room temperature) and 
To = 301 K we obtain a critical stress in the range ac = 
0.4-1.1 GPa, i.e. CTc = G/90-G/30. For comparison, we 
note that the shear yield strength of the metallic glass 
Vitreloy 1 is about 0.8 GPa = G/40 [3,[il, a value con- 
sistent with our estimate of the critical stress necessary 
for thermal runaway. 

The magnitude of the critical stress in these two esti- 
mates is remarkably similar considering the very different 
types of materials discussed and the enormous difference 
in time and length scales. However, as is evident from 
eq. ([6]), the quantities which govern the kinetic processes 
in these systems (in particular the scale h and the poorly 
constrained viscosity fi) do not appear in the expression 



for (Tc. It is not surprising, therefore, that the stress 
required to initiate spontaneous thermal runaway is rel- 
atively well constrained to about 1 GPa. 

These estimates and the concept of self-localizing ther- 
mal runaway demonstrate that our simple model is suffi- 
cient to explain failure below Frenkel's ultimate shear 
strength limit and strain localization at scales much 
larger than the lattice spacing. The fact that initiation of 
thermal runaway depends so weakly on the kinetic quan- 
tities gives confidence in the application of such models 
even to the very l arg e scales involved in continental de- 
formation 
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